
packages = c("brglm","ggcorrplot","interplot", "stargazer", "tidyverse",
             "here","haven","mice","modelsummary", "sandwich", "lmtest",
             "MASS", "ROCR", "generalhoslem", "brant", "bcROCsurface", "gridExtra",
             "erer", "margins", "interplot.medline")
#Warning! This will install R packages in your system if not yet installed. 
#load or install & load all packages 
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)



LEDs_data_f_mo <- read.csv(here::here("replication_data_missingobs_afs.csv"))
LEDs_data_f <- read.csv(here::here("replication_data_afs.csv"))


################################################################################# 
#                     Replications                                             #
################################################################################# 

#################### Summary Statistics ##################### 

stargazer(LEDs_data_f)
stargazer(LEDs_data_f_mo)

LEDs_data_f_mo$wdldownes <- as.factor(LEDs_data_f_mo$wdldownes)
LEDs_data_f$wdldownes <- as.factor(LEDs_data_f$wdldownes)

#Corr table for key variables
#Table 1
LEDs_data_corr_table <- LEDs_data_f %>% 
  ungroup() %>%
  mutate(wdldownes = as.numeric(wdldownes),
         wl = as.numeric(wl),
         logLER = as.numeric(logLER),
         wldownes = as.numeric(wldownes),) %>%
dplyr::select(wdldownes, wl,logLER,wldownes, relative_corr, initally, targally, v2x_libdem, v2x_polyarchy, politics , 
              v2x_ex_military, concap , capasst , qualrat , terrain , straterr , strat1 , strat2 , strat3 , strat4, 
              gdppc_growth, e_gdppc, relative_gdppc, relative_gdppc_growth, v2regsupgroups_6) %>%
  drop_na() 

corr <- cor(LEDs_data_corr_table)
colnames(corr) <- c("Lose/Draw/Win (D 2009)","Lose/Win (R&S 2002)", "LogLER","Lose/Win (D2009)", "Relative Corruption",
                    "Initiation", "Target",  "Liberal Democracy", "Poliarchy", "Polity", "Military Executive Dependence",
                    "Capabilities", "Alliance Capabilities", "Troop quality", "Terrain", "Strategy*Terrain",
                    "Strategy 1", "Strategy 2", "Strategy 3", "Strategy 4", 
                    "GDPPC Growth", "GDPPC", "Relative GDPPC", "Relative GDPPC Growth", "Ethnic Dependence" )
rownames(corr) <- c("Lose/Draw/Win (D 2009)","Lose/Win (R&S 2002)", "LogLER","Lose/Win (D2009)", "Relative Corruption",
                    "Initiation", "Target",  "Liberal Democracy", "Poliarchy", "Polity", "Military Executive Dependence",
                    "Capabilities", "Alliance Capabilities", "Troop quality", "Terrain", "Strategy*Terrain",
                    "Strategy 1", "Strategy 2", "Strategy 3", "Strategy 4", 
                    "GDPPC Growth", "GDPPC", "Relative GDPPC", "Relative GDPPC Growth", "Ethnic Dependence")
p.mat <- cor_pmat(LEDs_data_corr_table)
colnames(p.mat) <- c("Lose/Draw/Win (D 2009)","Lose/Win (R&S 2002)", "LogLER","Lose/Win (D2009)", "Relative Corruption",
                    "Initiation", "Target",  "Liberal Democracy", "Poliarchy", "Polity", "Military Executive Dependence",
                    "Capabilities", "Alliance Capabilities", "Troop quality", "Terrain", "Strategy*Terrain",
                    "Strategy 1", "Strategy 2", "Strategy 3", "Strategy 4", 
                    "GDPPC Growth", "GDPPC", "Relative GDPPC", "Relative GDPPC Growth", "Ethnic Dependence" )
rownames(p.mat) <- c("Lose/Draw/Win (D 2009)","Lose/Win (R&S 2002)", "LogLER","Lose/Win (D2009)", "Relative Corruption",
                    "Initiation", "Target",  "Liberal Democracy", "Poliarchy", "Polity", "Military Executive Dependence",
                    "Capabilities", "Alliance Capabilities", "Troop quality", "Terrain", "Strategy*Terrain",
                    "Strategy 1", "Strategy 2", "Strategy 3", "Strategy 4", 
                    "GDPPC Growth", "GDPPC", "Relative GDPPC", "Relative GDPPC Growth", "Ethnic Dependence")
ggcorrplot(corr, hc.order = F, p.mat = p.mat, sig.level = .05)
stargazer(cor(LEDs_data_corr_table))

####################Democracies at War (Reiter and Stam 2002)##################### 

iRSm1 <- glm(wl ~ polini + poltarg + init + concap + capasst + qualrat + 
               terrain + straterr + strat1 + strat2 + strat3 + strat4, 
             family = binomial(link = "probit"), data = LEDs_data_f)
summary(iRSm1)
iRSm1sr <- coeftest(iRSm1, 
                    vcov = vcovHC,
                    type = "HC1") #stata robust 

iRSm2 <- glm(wl ~ relative_corr + polini + poltarg + init + concap + capasst + qualrat + 
               terrain + straterr + strat1 + strat2 + strat3 + strat4,
             family = binomial(link = "probit"), data = LEDs_data_f)
summary(iRSm2)
iRSm2sr <- coeftest(iRSm2, 
                    vcov = vcovHC,
                    type = "HC1") #stata robust 

iRSm3 <- glm(wl ~ relative_corr_init + relative_corr_target + polini + poltarg + init + concap + capasst + qualrat + 
               terrain + straterr + strat1 + strat2 + strat3 + strat4, 
             family = binomial(link = "probit"), data = LEDs_data_f)
summary(iRSm3)
iRSm3sr <- coeftest(iRSm3, 
                    vcov = vcovHC,
                    type = "HC1") #stata robust 

#with LERD
iRSCLm1 <- lm(logLER ~ polini + poltarg + init + concap + capasst + qualrat + 
                terrain + straterr + strat1 + strat2 + strat3 + strat4, data = LEDs_data_f) 
summary(iRSCLm1)
iRSCLm1sr <- coeftest(iRSCLm1, 
                      vcov = vcovHC,
                      type = "HC1") #stata robust 

iRSCLm2 <- lm(logLER ~  relative_corr + polini + poltarg + init + concap + capasst + qualrat + 
                terrain + straterr + strat1 + strat2 + strat3 + strat4 , data = LEDs_data_f) 
summary(iRSCLm2)
iRSCLm2sr <- coeftest(iRSCLm2, 
                      vcov = vcovHC,
                      type = "HC1") #stata robust

iRSCLm3 <- lm(logLER ~ relative_corr_init+ relative_corr_target + 
                polini + poltarg + init + concap + capasst + qualrat + 
                terrain + straterr + strat1 + strat2 + strat3 + strat4, data = LEDs_data_f) 
summary(iRSCLm3)
iRSCLm3sr <- coeftest(iRSCLm3, 
                      vcov = vcovHC,
                      type = "HC1") #stata robust 

stargazer(iRSm1, iRSm2, iRSm3, iRSCLm1, iRSCLm2, iRSCLm3,
          se = list(iRSm1sr[,2], iRSm2sr[,2], iRSm3sr[,2],
                    iRSCLm1sr[,2], iRSCLm2sr[,2],iRSCLm3sr[,2]),
          covariate.labels = c("Relative Corruption", "Relative Corruption*Initiation", "Relative Corruption*Target",
                               "Politics*Initiation", "Politics*Target", "Initiation ",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c("Win/Lose", 
                            "Cochran and Long (2017) Battle Performance"))

          
##############How smart and tough are democracies? Downes (2009)##################


iDm1<- polr(wdldownes ~ pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f, method = c("probit"))
summary(iDm1)
iDm1sr <- coeftest(iDm1, 
                   vcov  = sandwich(x = iDm1, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict

iDm2<- polr(wdldownes ~ relative_corr + pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f, method = c("probit"))
summary(iDm2)
iDm2sr <- coeftest(iDm2, 
                   vcov = sandwich(x =iDm2, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict

iDm3<- polr(wdldownes ~ relative_corr_init + relative_corr_target
            + pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f, method = c("probit"))
summary(iDm3)
iDm3sr <- coeftest(iDm3, 
                   vcov = sandwich(x =iDm3, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict
#with LEDs

iDCLm1 <- lm(logLER ~ pol21+ initally+ targally+ 
               pol21initally+ pol21targally+ concap+ capasst+
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f) 
summary(iDCLm1)
iDCLm1sr <- coeftest(iDCLm1, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war)  #stata robust 

iDCLm2 <- lm(logLER ~ relative_corr + pol21+ initally+ targally+ 
               pol21initally+ pol21targally+ concap+ capasst+
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f) 
summary(iDCLm2)
iDCLm2sr <- coeftest(iDCLm2, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war)  #stata robust  

iDCLm3<- lm(logLER ~ relative_corr_init + relative_corr_target + 
              pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f)
summary(iDCLm3)
iDCLm3sr <- coeftest(iDCLm3, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war) #stata robust, slightly different because of PC rather than HC method. Clustered by conflict

#With Win Lose
iDWLm1 <- glm(wldownes ~ pol21+ initally+ targally+ 
               pol21initally+ pol21targally+ concap+ capasst+
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f, family = binomial(link = "probit")) 
summary(iDWLm1)
iDWLm1sr <- coeftest(iDWLm1, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war)  #stata robust 

iDWLm2 <- glm(wldownes ~ relative_corr + pol21+ initally+ targally+ 
               pol21initally+ pol21targally+ concap+ capasst+
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f, family = binomial(link = "probit")) 
summary(iDWLm2)
iDWLm2sr <- coeftest(iDWLm2, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war)  #stata robust  

iDWLm3<- glm(wldownes ~ relative_corr_init + relative_corr_target + 
              pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f, family = binomial(link = "probit"))
summary(iDWLm3)
iDWLm3sr <- coeftest(iDWLm3, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict


stargazer(iDm1, iDm2, iDm3, iDWLm1, iDWLm2,iDWLm3,iDCLm1, iDCLm2, iDCLm3,
          se = list(iDm1sr[,2], iDm2sr[,2], iDm3sr[,2],
                    iDWLm1sr[,2], iDWLm2sr[,2],iDWLm3sr[,2],
                    iDCLm1sr[,2], iDCLm2sr[,2], iDCLm3sr[,2]),
          covariate.labels = c("Relative Corruption", "Relative Corruption*Initiation", "Relative Corruption*Target",
                               "Polity", "Initiation", "Target",
                               "Polity*Initiation", "Polity*Target",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c("Win/Draw/Lose", 
                            'Win/Lose',
                            "Cochran and Long (2017) Battle Performance")
)

############Another skirmish over democracies (Reiter and Stam 2009)#############

iRS2m1 <- glm(wl ~ pol21targally + initally + pcini1old + pcini2old + concap +
                capasst + qualrat + terrain + straterr + strat1 + strat2 + strat3 + strat4, 
              family = binomial(link = "probit"), data = LEDs_data_f)
summary(iRS2m1)
iRS2m1sr <- coeftest(iRS2m1, 
                     vcov = vcovHC,
                     type = "HC1") #stata robust 


iRS2m2 <- glm(wl ~ relative_corr + pol21targally + initally + pcini1old + pcini2old + concap +
                capasst + qualrat + terrain + straterr + strat1 + strat2 + strat3 + strat4,
              family = binomial(link = "probit"), data = LEDs_data_f)
summary(iRS2m2)
iRS2m2sr <-coeftest(iRS2m2, 
                    vcov = vcovHC,
                    type = "HC1") #stata robust 

#with LERD
iRSCL2m1 <- lm(logLER ~ pol21targally + initally + pcini1old + pcini2old + concap +
                 capasst + qualrat + terrain + straterr + strat1 + strat2 + strat3 + strat4,
               data = LEDs_data_f) 
summary(iRSCL2m1)
iRSCL2m1sr <- coeftest(iRSCL2m1, 
                       vcov = vcovHC,
                       type = "HC1") #stata robust 

iRSCL2m2 <- lm(logLER ~  relative_corr + pol21targally + initally + pcini1old + pcini2old + concap +
                 capasst + qualrat + terrain + straterr + strat1 + strat2 + strat3 + strat4, data = LEDs_data_f) 
summary(iRSCL2m2)
iRSCL2m2sr <- coeftest(iRSCL2m2, 
                       vcov = vcovHC,
                       type = "HC1") #stata robust

stargazer(iRS2m1, iRS2m2, iRSCL2m1 ,iRSCL2m2,
          se = list(iRS2m1sr[,2], iRS2m2sr[,2],
                    iRSCL2m1sr[,2], iRSCL2m2sr[,2]),
          covariate.labels = c("Relative Corruption", "Politics*Target","Initiation",
                               "Politics (first curvilinear term)",
                               "Politics (second curvilinear term)",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c("Win/Lose", 
                            "Cochran and Long (2017) Battle Performance")
          
)

################################################################################# 
#                     Non-Imputed Data                                          #
################################################################################# 

#######Our own models########

m1_mo<- polr(wdldownes ~ relative_corr +
            initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
            qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
          data = LEDs_data_f_mo, method = c("probit"))
summary(m1_mo)
m1sr_mo<- coeftest(m1_mo, 
                vcov = sandwich(x =m1_mo, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_mo<- brglm(wl ~ relative_corr +
           initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f_mo,  family = binomial(link = "probit"), maxit = 200)
summary(m2_mo)
m2sr_mo<- coeftest(m2_mo, 
                vcov = sandwich(x =m2_mo, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_mo <- lm(logLER ~ relative_corr +
           + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f_mo) 
summary(m3_mo)
m3sr_mo<- coeftest(m3_mo, 
                vcov = vcovHC,
                type = "HC1",
                cluster = war) #stata robust 

m4_mo<- glm(wldownes ~ relative_corr +
           initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f_mo,  family = binomial(link = "probit"))
summary(m4_mo)
m4sr_mo<- coeftest(m4_mo, 
                vcov = sandwich(x = m4_mo, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_mo, m2_mo,m4_mo, m3_mo,
          se = list(m1sr_mo[,2], m2sr_mo[,2],m4sr_mo[,2], m3sr_mo[,2]),
          covariate.labels = c("Relative Corruption", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))


####################Democracies at War (Reiter and Stam 2002)##################### 

iRSm1_mo<- brglm(wl ~ polini + poltarg + init + concap + capasst + qualrat + 
                 terrain + straterr + strat1 + strat2 + strat3 + strat4, 
               family = binomial(link = "probit"), data = LEDs_data_f_mo, maxit= 200)
summary(iRSm1_mo)
iRSm1sr_mo <- coeftest(iRSm1_mo, 
                    vcov = vcovHC,
                    type = "HC1") #stata robust 

iRSm2_mo<- brglm(wl ~ relative_corr + polini + poltarg + init + concap + capasst + qualrat + 
                 terrain + straterr + strat1 + strat2 + strat3 + strat4,
               family = binomial(link = "probit"), data = LEDs_data_f_mo )
summary(iRSm2_mo)
iRSm2sr_mo <- coeftest(iRSm2_mo, 
                    vcov = vcovHC,
                    type = "HC1") #stata robust 

iRSm3_mo<- brglm(wl ~ relative_corr_init + relative_corr_target + polini + poltarg + init + concap + capasst + qualrat + 
                 terrain + straterr + strat1 + strat2 + strat3 + strat4, 
               family = binomial(link = "probit"), data = LEDs_data_f_mo, maxit=500 )
summary(iRSm3_mo)
iRSm3sr_mo <- coeftest(iRSm3_mo, 
                    vcov = vcovHC,
                    type = "HC1") #stata robust 

#with LERD
iRSCLm1_mo<- lm(logLER ~ polini + poltarg + init + concap + capasst + qualrat + 
                  terrain + straterr + strat1 + strat2 + strat3 + strat4, data = LEDs_data_f_mo) 
summary(iRSCLm1_mo)
iRSCLm1sr_mo <- coeftest(iRSCLm1_mo, 
                      vcov = vcovHC,
                      type = "HC1") #stata robust 

iRSCLm2_mo<- lm(logLER ~  relative_corr + polini + poltarg + init + concap + capasst + qualrat + 
                  terrain + straterr + strat1 + strat2 + strat3 + strat4 , data = LEDs_data_f_mo) 
summary(iRSCLm2_mo)
iRSCLm2sr_mo <- coeftest(iRSCLm2_mo, 
                      vcov = vcovHC,
                      type = "HC1") #stata robust

iRSCLm3_mo<- lm(logLER ~ relative_corr_init+ relative_corr_target + 
                  polini + poltarg + init + concap + capasst + qualrat + 
                  terrain + straterr + strat1 + strat2 + strat3 + strat4, data = LEDs_data_f_mo) 
summary(iRSCLm3_mo)
iRSCLm3sr_mo <- coeftest(iRSCLm3_mo, 
                      vcov = vcovHC,
                      type = "HC1") #stata robust 

stargazer(iRSm1_mo, iRSm2_mo, iRSm3_mo, iRSCLm1_mo, iRSCLm2_mo, iRSCLm3_mo,
          se = list(iRSm1sr_mo[,2], iRSm2sr_mo[,2], iRSm3sr_mo[,2],
                    iRSCLm1sr_mo[,2], iRSCLm2sr_mo[,2],iRSCLm3sr_mo[,2]),
          covariate.labels = c("Relative Corruption", "Relative Corruption*Initiation", "Relative Corruption*Target",
                               "Politics*Initiation", "Politics*Target",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c("Win/Lose", 
                            "Cochran and Long (2017) Battle Performance"))

##############How smart and tough are democracies? Downes (2009)##################

iDm1_mo<- polr(wdldownes ~ pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f_mo, method = c("probit"))
summary(iDm1_mo)
iDm1sr_mo <- coeftest(iDm1_mo, 
                   vcov  = sandwich(x = iDm1_mo, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict

iDm2_mo<- polr(as.factor(wdldownes) ~ relative_corr + pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f_mo, method = c("probit"))
summary(iDm2_mo)
iDm2sr_mo <- coeftest(iDm2_mo, 
                   vcov = sandwich(x =iDm2_mo, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict

iDm3_mo<- polr(wdldownes ~ relative_corr_init + relative_corr_target
               + pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f_mo, method = c("probit"))
summary(iDm3_mo)
iDm3sr_mo <- coeftest(iDm3_mo, 
                   vcov = sandwich(x =iDm3_mo, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict
#with LEDs

iDCLm1_mo<- lm(logLER ~ pol21+ initally+ targally+ 
                 pol21initally+ pol21targally+ concap+ capasst+
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f_mo)
summary(iDCLm1_mo)
iDCLm1sr_mo <- coeftest(iDCLm1_mo, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war)  #stata robust 

iDCLm2_mo<- lm(logLER ~ relative_corr + pol21+ initally+ targally+ 
                 pol21initally+ pol21targally+ concap+ capasst+
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f_mo) 
summary(iDCLm2_mo)
iDCLm2sr_mo <- coeftest(iDCLm2_mo, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war)  #stata robust  

iDCLm3_mo<- lm(logLER ~ relative_corr_init + relative_corr_target + 
                 pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f_mo)
summary(iDCLm3_mo)
iDCLm3sr_mo <- coeftest(iDCLm3_mo, 
                     vcov = sandwich(x = iDCLm3_mo, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict


#With Win Lose
iDWLm1_mo <- glm(wldownes ~ pol21+ initally+ targally+ 
                pol21initally+ pol21targally+ concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f_mo, family = binomial(link = "probit")) 
summary(iDWLm1_mo)
iDWLm1sr_mo <- coeftest(iDWLm1_mo, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war)  #stata robust 

iDWLm2_mo <- glm(wldownes ~ relative_corr + pol21+ initally+ targally+ 
                pol21initally+ pol21targally+ concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f_mo, family = binomial(link = "probit")) 
summary(iDWLm2_mo)
iDWLm2sr_mo <- coeftest(iDWLm2_mo, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war)  #stata robust  

iDWLm3_mo<- glm(wldownes ~ relative_corr_init + relative_corr_target + 
               pol21+ initally+ targally+ pol21initally+ pol21targally+ concap+ capasst+
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f_mo, family = binomial(link = "probit"))
summary(iDWLm3_mo)
iDWLm3sr_mo <- coeftest(iDWLm3_mo, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war) #stata robust, slightly different because of of PC rather than HC method. Clustered by conflict

stargazer(iDm1_mo, iDm2_mo, iDm3_mo,iDWLm1_mo, iDWLm2_mo, iDWLm3_mo, iDCLm1_mo, iDCLm2_mo, iDCLm3_mo,
          se = list(iDm1sr_mo[,2], iDm2sr_mo[,2], iDm3sr_mo[,2],
                    iDWLm1sr_mo[,2],iDWLm2sr_mo[,2],iDWLm3sr_mo[,2],
                    iDCLm1sr_mo[,2], iDCLm2sr_mo[,2], iDCLm3sr_mo[,2]),
          covariate.labels = c("Relative Corruption", "Relative Corruption*Initiation", "Relative Corruption*Target",
                               "Polity", "Initiation", "Target",
                               "Polity*Initiation", "Polity*Target",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c("Win/Draw/Lose", 
                            'Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))
          
############Another skirmish over democracies (Reiter and Stam 2009)#############

iRS2m1_mo<- brglm(wl ~ pol21targally + initally + pcini1old + pcini2old + concap +
                  capasst + qualrat + terrain + straterr + strat1 + strat2 + strat3 + strat4, 
                family = binomial(link = "probit"), data = LEDs_data_f_mo)
summary(iRS2m1_mo)
iRS2m1sr_mo <- coeftest(iRS2m1_mo, 
                     vcov = vcovHC,
                     type = "HC1") #stata robust 


iRS2m2_mo<- brglm(wl ~ relative_corr + pol21targally + initally + pcini1old + pcini2old + concap +
                  capasst + qualrat + terrain + straterr + strat1 + strat2 + strat3 + strat4,
                family = binomial(link = "probit"), data = LEDs_data_f_mo)
summary(iRS2m2_mo)
iRS2m2sr_mo <-coeftest(iRS2m2_mo, 
                    vcov = vcovHC,
                    type = "HC1") #stata robust 

#with LERD
iRSCL2m1_mo<- lm(logLER ~ pol21targally + initally + pcini1old + pcini2old + concap +
                   capasst + qualrat + terrain + straterr + strat1 + strat2 + strat3 + strat4,
                 data = LEDs_data_f_mo) 
summary(iRSCL2m1_mo)
iRSCL2m1sr_mo <- coeftest(iRSCL2m1_mo, 
                       vcov = vcovHC,
                       type = "HC1") #stata robust 

iRSCL2m2_mo<- lm(logLER ~  relative_corr + pol21targally + initally + pcini1old + pcini2old + concap +
                   capasst + qualrat + terrain + straterr + strat1 + strat2 + strat3 + strat4, data = LEDs_data_f_mo) 
summary(iRSCL2m2_mo)
iRSCL2m2sr_mo <- coeftest(iRSCL2m2_mo, 
                       vcov = vcovHC,
                       type = "HC1") #stata robust

stargazer(iRS2m1_mo, iRS2m2_mo, iRSCL2m1_mo,iRSCL2m2_mo,
          se = list(iRS2m1sr_mo[,2], iRS2m2sr_mo[,2],
                    iRSCL2m1sr_mo[,2], iRSCL2m2sr_mo[,2]),
          covariate.labels = c("Relative Corruption", "Politics*Target","Initiation",
                               "Politics (first curvilinear term)",
                               "Politics (second curvilinear term)",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c("Win/Lose", 
                            "Cochran and Long (2017) Battle Performance")
          
)


################################################################################# 
#                     Different Model Specifications                            #
################################################################################# 

#######No Control Variables########

m1_nc<- polr(as.factor(wdldownes) ~ relative_corr, 
             data = LEDs_data_f, method = c("probit"))
summary(m1_nc)
m1sr_nc<- coeftest(m1_nc, 
                   vcov = sandwich(x =m1_nc, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_nc<- glm(wl ~ relative_corr,
            data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_nc)
m2sr_nc<- coeftest(m2_nc, 
                   vcov = sandwich(x =m2_nc, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_nc <- lm(logLER ~ relative_corr, 
            data = LEDs_data_f) 
summary(m3_nc)
m3sr_nc<- coeftest(m3_nc, 
                   vcov = vcovHC,
                   type = "HC1",
                   cluster = war) #stata robust 

m4_nc<- glm(wldownes ~ relative_corr, 
            data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_nc)
m4sr_nc<- coeftest(m4_nc,
                   vcov = sandwich(x = m4_nc, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_nc, m2_nc,m4_nc, m3_nc,
          se = list(m1sr_nc[,2], m2sr_nc[,2],m4sr_nc[,2], m3sr_nc[,2]),
          covariate.labels = c("Relative Corruption"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

#######Country Level Control Variables########

m1_clc<- polr(as.factor(wdldownes) ~ relative_corr +
                v2x_libdem+  concap + capasst + qualrat, 
              data = LEDs_data_f, method = c("probit"))
summary(m1_clc)
m1sr_clc<- coeftest(m1_clc, 
                    vcov = sandwich(x =m1_clc, adjust = TRUE),
                    cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_clc<- glm(wl ~ relative_corr +
               v2x_libdem+  concap + capasst + qualrat, 
             data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_clc)
m2sr_clc<- coeftest(m2_clc, 
                    vcov = sandwich(x =m2_clc, adjust = TRUE),
                    cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_clc <- lm(logLER ~ relative_corr +
               v2x_libdem+  concap + capasst + qualrat, 
             data = LEDs_data_f) 
summary(m3_clc)
m3sr_clc<- coeftest(m3_clc, 
                    vcov = vcovHC,
                    type = "HC1",
                    cluster = war) #stata robust 

m4_clc<- glm(wldownes ~ relative_corr +
               v2x_libdem+  concap + capasst + qualrat, 
             data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_clc)
m4sr_clc<- coeftest(m4_clc, 
                    vcov = sandwich(x = m4_clc, adjust = TRUE),
                    cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_clc, m2_clc,m4_clc, m3_clc,
          se = list(m1sr_clc[,2], m2sr_clc[,2],m4sr_clc[,2], m3sr_clc[,2]),
          covariate.labels = c("Relative Corruption", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))


#######Extra Control Variables########

####### GDPPC ########
m1_gdp<- polr(as.factor(wdldownes) ~ relative_corr +  e_gdppc + 
               initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f, method = c("probit"))
summary(m1_gdp)
m1sr_gdp<- coeftest(m1_gdp, 
                   vcov = sandwich(x =m1_gdp, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_gdp<- glm(wl ~ relative_corr + e_gdppc  + 
              initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_gdp)
m2sr_gdp<- coeftest(m2_gdp, 
                   vcov = sandwich(x =m2_gdp, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_gdp <- lm(logLER ~ relative_corr + e_gdppc  + 
              + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f) 
summary(m3_gdp)
m3sr_gdp<- coeftest(m3_gdp, 
                   vcov = vcovHC,
                   type = "HC1",
                   cluster = war) #stata robust 

m4_gdp<- glm(wldownes ~ relative_corr + e_gdppc  + 
              initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_gdp)
m4sr_gdp<- coeftest(m4_gdp, 
                   vcov = sandwich(x = m4_gdp, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_gdp, m2_gdp,m4_gdp, m3_gdp,
          se = list(m1sr_gdp[,2], m2sr_gdp[,2],m4sr_gdp[,2], m3sr_gdp[,2]),
          covariate.labels = c("Relative Corruption", "GDPPC", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

####### GDPPC Growth ########
m1_gdpg<- polr(as.factor(wdldownes) ~ relative_corr +  gdppc_growth + 
                initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f, method = c("probit"))
summary(m1_gdpg)
m1sr_gdpg<- coeftest(m1_gdpg, 
                    vcov = sandwich(x =m1_gdpg, adjust = TRUE),
                    cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_gdpg<- glm(wl ~ relative_corr + gdppc_growth  + 
               initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_gdpg)
m2sr_gdpg<- coeftest(m2_gdpg, 
                    vcov = sandwich(x =m2_gdpg, adjust = TRUE),
                    cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_gdpg <- lm(logLER ~ relative_corr + gdppc_growth  + 
               + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f) 
summary(m3_gdpg)
m3sr_gdpg<- coeftest(m3_gdpg, 
                    vcov = vcovHC,
                    type = "HC1",
                    cluster = war) #stata robust 

m4_gdpg<- glm(wldownes ~ relative_corr + gdppc_growth  + 
               initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_gdpg)
m4sr_gdpg<- coeftest(m4_gdpg, 
                    vcov = sandwich(x = m4_gdpg, adjust = TRUE),
                    cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_gdpg, m2_gdpg,m4_gdpg, m3_gdpg,
          se = list(m1sr_gdpg[,2], m2sr_gdpg[,2],m4sr_gdpg[,2], m3sr_gdpg[,2]),
          covariate.labels = c("Relative Corruption", "GDPPC Growth", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

####### Relative GDPPC #######
m1_gdpg_re<- polr(as.factor(wdldownes) ~ relative_corr +  relative_gdppc + 
                 initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f, method = c("probit"))
summary(m1_gdpg_re)
m1sr_gdpg_re<- coeftest(m1_gdpg_re, 
                     vcov = sandwich(x =m1_gdpg_re, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_gdpg_re<- glm(wl ~ relative_corr + relative_gdppc  + 
                initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_gdpg_re)
m2sr_gdpg_re<- coeftest(m2_gdpg_re, 
                     vcov = sandwich(x =m2_gdpg_re, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_gdpg_re <- lm(logLER ~ relative_corr + relative_gdppc  + 
                + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f) 
summary(m3_gdpg_re)
m3sr_gdpg_re<- coeftest(m3_gdpg_re, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war) #stata robust 

m4_gdpg_re<- glm(wldownes ~ relative_corr + relative_gdppc  + 
                initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_gdpg_re)
m4sr_gdpg_re<- coeftest(m4_gdpg_re, 
                     vcov = sandwich(x = m4_gdpg_re, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_gdpg_re, m2_gdpg_re,m4_gdpg_re, m3_gdpg_re,
          se = list(m1sr_gdpg_re[,2], m2sr_gdpg_re[,2],m4sr_gdpg_re[,2], m3sr_gdpg_re[,2]),
          covariate.labels = c("Relative Corruption", "Relative GDPPC", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))


####### Relative GDP Growth ####
m1_gdpgg_re<- polr(as.factor(wdldownes) ~ relative_corr +  relative_gdppc_growth + 
                    initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                    qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                  data = LEDs_data_f, method = c("probit"))
summary(m1_gdpgg_re)
m1sr_gdpgg_re<- coeftest(m1_gdpgg_re, 
                        vcov = sandwich(x =m1_gdpgg_re, adjust = TRUE),
                        cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_gdpgg_re<- glm(wl ~ relative_corr + relative_gdppc_growth  + 
                   initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                   qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                 data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_gdpgg_re)
m2sr_gdpgg_re<- coeftest(m2_gdpgg_re, 
                        vcov = sandwich(x =m2_gdpgg_re, adjust = TRUE),
                        cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_gdpgg_re <- lm(logLER ~ relative_corr + relative_gdppc_growth  + 
                   + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                   qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                 data = LEDs_data_f) 
summary(m3_gdpgg_re)
m3sr_gdpgg_re<- coeftest(m3_gdpgg_re, 
                        vcov = vcovHC,
                        type = "HC1",
                        cluster = war) #stata robust 

m4_gdpgg_re<- glm(wldownes ~ relative_corr + relative_gdppc_growth  + 
                   initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                   qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                 data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_gdpgg_re)
m4sr_gdpgg_re<- coeftest(m4_gdpgg_re, 
                        vcov = sandwich(x = m4_gdpgg_re, adjust = TRUE),
                        cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_gdpgg_re, m2_gdpgg_re,m4_gdpgg_re, m3_gdpgg_re,
          se = list(m1sr_gdpgg_re[,2], m2sr_gdpgg_re[,2],m4sr_gdpgg_re[,2], m3sr_gdpgg_re[,2]),
          covariate.labels = c("Relative Corruption", "Relative GDPPC Growth", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

######## Relative GDPPC Growth + Ethnic Dependency #######

m1_av<- polr(as.factor(wdldownes) ~ relative_corr +  relative_gdppc_growth + v2regsupgroups_6 + 
               initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_f, method = c("probit"))
summary(m1_av)
m1sr_av<- coeftest(m1_av, 
                   vcov = sandwich(x =m1_av, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_av<- glm(wl ~ relative_corr + relative_gdppc_growth + v2regsupgroups_6 + 
              initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_av)
m2sr_av<- coeftest(m2_av, 
                   vcov = sandwich(x =m2_av, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_av <- lm(logLER ~ relative_corr + relative_gdppc_growth + v2regsupgroups_6 + 
              + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f) 
summary(m3_av)
m3sr_av<- coeftest(m3_av, 
                   vcov = vcovHC,
                   type = "HC1",
                   cluster = war) #stata robust 

m4_av<- glm(wldownes ~ relative_corr + relative_gdppc_growth + v2regsupgroups_6 + 
              initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_av)
m4sr_av<- coeftest(m4_av, 
                   vcov = sandwich(x = m4_av, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_av, m2_av,m4_av, m3_av,
          se = list(m1sr_av[,2], m2sr_av[,2],m4sr_av[,2], m3sr_av[,2]),
          covariate.labels = c("Relative Corruption", "Relative GDPPC Growth", "Ethnic Dependence", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

#######No Military Executive Dependency#######
m1_noed<- polr(as.factor(wdldownes) ~ relative_corr +
            initally + targally + v2x_libdem + concap + capasst +
            qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
          data = LEDs_data_f, method = c("probit"))
summary(m1_noed)
ocME(w = m1_noed)
m1sr_noed<- coeftest(m1_noed, 
                vcov = sandwich(x =m1_noed, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_noed<- glm(wl ~ relative_corr +
           initally + targally + v2x_libdem + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_noed)
m2sr_noed<- coeftest(m2_noed, 
                vcov = sandwich(x =m2_noed, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m3_noed<- glm(wldownes ~ relative_corr +
           initally + targally + v2x_libdem + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m3_noed)
m3sr_noed<- coeftest(m3_noed, 
                vcov = sandwich(x = m3_noed, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m4_noed <- lm(logLER ~ relative_corr +
           + initally +targally + v2x_libdem + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f) 
summary(m4_noed)
m4sr_noed<- coeftest(m4_noed, 
                vcov = vcovHC,
                type = "HC1",
                cluster = war) #stata robust 

stargazer(m1_noed, m2_noed, m3_noed, m4_noed,
          se = list(m1sr_noed[,2], m2sr_noed[,2],m3sr_noed[,2], m4sr_noed[,2]),
          covariate.labels = c("Corruption", "Initiation", "Target", 
                               "Liberal Democracy",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

####### Polyarchy ########

m1_poly<- polr(as.factor(wdldownes) ~ relative_corr +
            initally + targally + v2x_polyarchy+ v2x_ex_military + concap + capasst +
            qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
          data = LEDs_data_f, method = c("probit"))
summary(m1_poly)
ocME(w = m1_poly)
m1sr_poly<- coeftest(m1_poly, 
                vcov = sandwich(x =m1_poly, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_poly<- glm(wl ~ relative_corr +
           initally + targally + v2x_polyarchy+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_poly)
m2sr_poly<- coeftest(m2_poly, 
                vcov = sandwich(x =m2_poly, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m3_poly<- glm(wldownes ~ relative_corr +
           initally + targally + v2x_polyarchy+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m3_poly)
m3sr_poly<- coeftest(m3_poly, 
                vcov = sandwich(x = m3_poly, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m4_poly <- lm(logLER ~ relative_corr +
           + initally +targally + v2x_polyarchy+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f) 
summary(m4_poly)
m4sr_poly<- coeftest(m4_poly, 
                vcov = vcovHC,
                type = "HC1",
                cluster = war) #stata robust 

stargazer(m1_poly, m2_poly, m3_poly, m4_poly,
          se = list(m1sr_poly[,2], m2sr_poly[,2],m3sr_poly[,2], m4sr_poly[,2]),
          covariate.labels = c("Corruption", "Initiation", "Target", 
                               "Polyarchy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

####### Polyarchy No Mil Dep ########

m1_polynmd<- polr(as.factor(wdldownes) ~ relative_corr +
                 initally + targally + v2x_polyarchy + concap + capasst +
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f, method = c("probit"))
summary(m1_polynmd)
ocME(w = m1_polynmd)
m1sr_polynmd<- coeftest(m1_polynmd, 
                     vcov = sandwich(x =m1_polynmd, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_polynmd<- glm(wl ~ relative_corr +
                initally + targally + v2x_polyarchy + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_polynmd)
m2sr_polynmd<- coeftest(m2_polynmd, 
                     vcov = sandwich(x =m2_polynmd, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

m3_polynmd<- glm(wldownes ~ relative_corr +
                initally + targally + v2x_polyarchy + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m3_polynmd)
m3sr_polynmd<- coeftest(m3_polynmd, 
                     vcov = sandwich(x = m3_polynmd, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m4_polynmd <- lm(logLER ~ relative_corr +
                + initally +targally + v2x_polyarchy + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f) 
summary(m4_polynmd)
m4sr_polynmd<- coeftest(m4_polynmd, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war) #stata robust 

stargazer(m1_polynmd, m2_polynmd, m3_polynmd, m4_polynmd,
          se = list(m1sr_polynmd[,2], m2sr_polynmd[,2],m3sr_polynmd[,2], m4sr_polynmd[,2]),
          covariate.labels = c("Corruption", "Initiation", "Target", 
                               "Polyarchy",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

####### Absolute Level of Corruption ########

m1_corr<- polr(as.factor(wdldownes) ~ v2x_corr +
                 initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f, method = c("probit"))
summary(m1_corr)
m1sr_corr<- coeftest(m1_corr, 
                     vcov = sandwich(x =m1_corr, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_corr<- glm(wl ~ v2x_corr +
                initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_corr)
m2sr_corr<- coeftest(m2_corr, 
                     vcov = sandwich(x =m2_corr, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_corr <- lm(logLER ~ v2x_corr +
                + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f) 
summary(m3_corr)
m3sr_corr<- coeftest(m3_corr, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war) #stata robust 

m4_corr<- glm(wldownes ~ v2x_corr +
                initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_corr)
m4sr_corr<- coeftest(m4_corr, 
                     vcov = sandwich(x = m4_corr, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_corr, m2_corr,m4_corr, m3_corr,
          se = list(m1sr_corr[,2], m2sr_corr[,2],m4sr_corr[,2], m3sr_corr[,2]),
          covariate.labels = c("Corruption", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

####### Average Oponent Level of Corruption ########

m1_avgcorr<- polr(as.factor(wdldownes) ~ opp_corr +
                    initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                    qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                  data = LEDs_data_f, method = c("probit"))
summary(m1_avgcorr)
m1sr_avgcorr<- coeftest(m1_avgcorr, 
                        vcov = sandwich(x =m1_avgcorr, adjust = TRUE),
                        cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_avgcorr<- glm(wl ~ opp_corr +
                   initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                   qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                 data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2_avgcorr)
m2sr_avgcorr<- coeftest(m2_avgcorr, 
                        vcov = sandwich(x =m2_avgcorr, adjust = TRUE),
                        cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_avgcorr <- lm(logLER ~ opp_corr +
                   + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                   qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                 data = LEDs_data_f) 
summary(m3_avgcorr)
m3sr_avgcorr<- coeftest(m3_avgcorr, 
                        vcov = vcovHC,
                        type = "HC1",
                        cluster = war) #stata robust 

m4_avgcorr<- glm(wldownes ~ opp_corr +
                   initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                   qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                 data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4_avgcorr)
m4sr_avgcorr<- coeftest(m4_avgcorr, 
                        vcov = sandwich(x = m4_avgcorr, adjust = TRUE),
                        cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_avgcorr, m2_avgcorr,m4_avgcorr, m3_avgcorr,
          se = list(m1sr_avgcorr[,2], m2sr_avgcorr[,2],m4sr_avgcorr[,2], m3sr_avgcorr[,2]),
          covariate.labels = c("Average Oponent Corruption", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

####################Main Belligerents####################

LEDs_data_mainbel <- LEDs_data_f %>% filter(initally == 1 | targally ==1 )


m1_mb<- polr(as.factor(wdldownes) ~ relative_corr +
               initally  + v2x_libdem+ v2x_ex_military + concap + capasst +
               qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
             data = LEDs_data_mainbel, method = c("probit"))
summary(m1_mb)
m1sr_mb<- coeftest(m1_mb, 
                   vcov = sandwich(x =m1_mb, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_mb<- glm(wl ~ relative_corr +
              initally  + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_mainbel,  family = binomial(link = "probit"))
summary(m2_mb)
m2sr_mb<- coeftest(m2_mb, 
                   vcov = sandwich(x =m2_mb, adjust = TRUE),
                   cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3_mb <- lm(logLER ~ relative_corr +
              + initally + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_mainbel) 
summary(m3_mb)
m3sr_mb<- coeftest(m3_mb, 
                   vcov = vcovHC,
                   type = "HC1",
                   cluster = war) #stata robust 

m4_mb<- glm(wldownes ~ relative_corr +
              initally + v2x_libdem+ v2x_ex_military + concap+ capasst+
              qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
            data = LEDs_data_mainbel,  family = binomial(link = "probit"))
summary(m4_mb)
m4sr_mb<- coeftest(m4_mb, 
                vcov = sandwich(x = m4_mb, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC


stargazer(m1_mb, m2_mb,m4_mb, m3_mb,
          se = list(m1sr_mb[,2], m2sr_mb[,2],m4sr_mb[,2], m3sr_mb[,2]),
          covariate.labels = c("Corruption", "Initiation",
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))

########### Bilateral Conflict ############
bi_lateral_conflict_data <- LEDs_data_f %>% group_by(RSwarnumber) %>% mutate(n_combatents = n()) %>% 
  filter(n_combatents == 2)

m1_bc<- polr(as.factor(wdldownes) ~ relative_corr +
            initally  + v2x_libdem+ v2x_ex_military + concap  +
            qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
          data = bi_lateral_conflict_data, method = c("probit"))
summary(m1_bc)
m1sr_bc<- coeftest(m1_bc, 
                vcov = sandwich(x =m1_bc, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m2_bc<- glm(wl ~ relative_corr +
           initally + v2x_libdem+ v2x_ex_military + concap+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = bi_lateral_conflict_data,  family = binomial(link = "probit"))
summary(m2_bc)
m2sr_bc<- coeftest(m2_bc, 
                vcov = sandwich(x =m2_bc, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m3_bc<- glm(wldownes ~ relative_corr +
           initally + v2x_libdem+ v2x_ex_military + concap+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = bi_lateral_conflict_data,  family = binomial(link = "probit"))
summary(m3_bc)
m3sr_bc<- coeftest(m3_bc, 
                vcov = sandwich(x = m3_bc, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m4_bc <- lm(logLER ~ relative_corr +
           + initally + v2x_libdem+ v2x_ex_military+ concap+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = bi_lateral_conflict_data) 
summary(m4_bc)
m4sr_bc<- coeftest(m4_bc, 
                vcov = vcovHC,
                type = "HC1",
                cluster = war) #stata robust 

stargazer(m1_bc, m2_bc, m3_bc, m4_bc,
          se = list(m1sr_bc[,2], m2sr_bc[,2],m3sr_bc[,2], m4sr_bc[,2]),
          covariate.labels = c("Corruption", "Initiation",
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))


################################################################################# 
#                         In Sample Model Fit Comparison                        #
################################################################################# 
########Re Run Main Models for Comparison#######
m1<- polr(wdldownes ~ relative_corr +
            initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
            qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
          data = LEDs_data_f, method = c("probit"))
summary(m1)
m1sr<- coeftest(m1, 
                vcov = sandwich(x =m1, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m2<- glm(wl ~ relative_corr +
           initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2)
m2sr<- coeftest(m2, 
                vcov = sandwich(x =m2, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m3 <- lm(logLER ~ relative_corr +
           + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f) 
summary(m3)
m3sr<- coeftest(m3, 
                vcov = vcovHC,
                type = "HC1",
                cluster = war) #stata robust 

m4<- glm(wldownes ~ relative_corr +
           initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m4)
m4sr<- coeftest(m4, 
                vcov = sandwich(x = m4, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC



################################################################################# 
#                                 Time Effects                                 #
################################################################################# 

LEDs_data_f$periods50 <- ifelse(LEDs_data_f$year <= 1873, 1,
                                ifelse(LEDs_data_f$year >1873 & LEDs_data_f$year <= 1923,2,
                                       ifelse(LEDs_data_f$year > 1923 & LEDs_data_f$year <= 1973,3,
                                              ifelse(LEDs_data_f$year > 1973, 4, NA))))
LEDs_data_f$periods1900 <- ifelse(LEDs_data_f$year <= 1900, 1,2)


#######Our own models 50 Year Periods########

im1_p50<- polr(as.factor(wdldownes) ~ relative_corr*periods50 + 
                 initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                 qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
               data = LEDs_data_f, method = c("probit"))
summary(im1_p50)
im1sr_p50<- coeftest(im1_p50, 
                     vcov = sandwich(x = im1_p50, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

#grab relative_corr from interplot.medline for compare extremes heuristic
var1 = "relative_corr"
m = im1_p50
var2 = "periods50"
coef1 <- coef(m)[var1]
coef3 <- coef(m)[ paste0(var1,":", var2)]
med <- median(eval(parse(text = paste0("im1_p50$model$", var2))))
hline <- coef1 + med*coef3

p50_1 <- interplot(im1_p50, var1 = "relative_corr", var2 = "periods50", hist = T, point = T)+
  labs(y="Corruption", x="Year", title = 'Downes (2009) Win/Draw/Lose') +
  ggplot2::geom_hline(yintercept=hline, linetype="dashed") +
  ggplot2::geom_hline(yintercept=0) +
  theme_bw()

im2_p50<- glm(wl ~ relative_corr*periods50 +
                initally + targally + v2x_libdem+ v2x_ex_military +concap + capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f,  family = binomial(link = "probit"))
summary(im2_p50)
im2sr_p50<- coeftest(im2_p50, 
                     vcov = sandwich(x = im2_p50, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC

#grab relative_corr from interplot.medline for compare extremes heuristic
var1 = "relative_corr"
m = im2_p50
var2 = "periods50"
coef1 <- coef(m)[var1]
coef3 <- coef(m)[ paste0(var1,":", var2)]
med <- median(eval(parse(text = paste0("im2_p50$model$", var2))))
hline <- coef1 + med*coef3

p50_2 <- interplot(im2_p50, var1 = "relative_corr", var2 = "periods50", hist = T, point = T)+
  labs(y="Corruption", x="Year", title = 'Reiter and Stam (2002) Win/Lose') +
  ggplot2::geom_hline(yintercept=hline, linetype="dashed") +
  ggplot2::geom_hline(yintercept=0) +
  theme_bw()


im3_p50<- glm(wldownes ~ relative_corr*periods50 +
                initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f,  family = binomial(link = "probit"))
summary(im3_p50)

im3sr_p50<- coeftest(im3_p50, 
                     vcov = sandwich(x = im3_p50, adjust = TRUE),
                     cluster = war) #stata robust, slightly different because of of PC rather than HC


#grab relative_corr from interplot.medline for compare extremes heuristic
var1 = "relative_corr"
m = im3_p50
var2 = "periods50"
coef1 <- coef(m)[var1]
coef3 <- coef(m)[ paste0(var1,":", var2)]
med <- median(eval(parse(text = paste0("im3_p50$model$", var2))))
hline <- coef1 + med*coef3

p50_3 <-interplot(im3_p50, var1 = "relative_corr", var2 = "periods50", hist = T, point = T)+
  labs(y="Corruption", x="Year", title =  "Downes (2009) Win/Lose") +
  ggplot2::geom_hline(yintercept=hline, linetype="dashed") +
  ggplot2::geom_hline(yintercept=0) +
  theme_bw()

#with LEDs

im4_p50 <- lm(logLER ~ relative_corr*periods50 +
                + initally +targally + v2x_libdem+ v2x_ex_military +concap+ capasst+
                qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
              data = LEDs_data_f) 
summary(im4_p50)
im4sr_p50<- coeftest(im4_p50, 
                     vcov = vcovHC,
                     type = "HC1",
                     cluster = war) #stata robust 

#grab relative_corr from interplot.medline for compare extremes heuristic
var1 = "relative_corr"
m = im4_p50
var2 = "periods50"
coef1 <- coef(m)[var1]
coef3 <- coef(m)[ paste0(var1,":", var2)]
med <- median(eval(parse(text = paste0("im4_p50$model$", var2))))
hline <- coef1 + med*coef3

p50_4 <- interplot(im4_p50, var1 = "relative_corr", var2 = "periods50", hist = T, point = T)+
  labs(y="Corruption", x="Year", title = "Cochran and Long (2017) Battle Performance") +
  ggplot2::geom_hline(yintercept=hline, linetype="dashed") +
  ggplot2::geom_hline(yintercept=0) +
  theme_bw()

stargazer(im1_p50, im2_p50,im3_p50, im4_p50,
          se = list(im1sr_p50[,2], im2sr_p50[,2], im3sr_p50[,2], im4sr_p50[,2]),
          covariate.labels = c("Relative Corruption", "Period(50 years)", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4", "Corruption*Period(50 years)"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose', 
                            "Downes (2009) Win/Lose",
                            "Cochran and Long (2017) Battle Performance"
          ))

grid.arrange(p50_1, p50_2, p50_3, p50_4, nrow = 2, bottom = "Model")

#######Our own models Post-1900########

im1_p1900<- polr(as.factor(wdldownes) ~ relative_corr*periods1900 + 
                   initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                   qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                 data = LEDs_data_f, method = c("probit"))
summary(im1_p1900)
im1sr_p1900<- coeftest(im1_p1900, 
                       vcov = sandwich(x = im1_p1900, adjust = TRUE),
                       cluster = war) #stata robust, slightly different because of of PC rather than HC


#grab relative_corr from interplot.medline for compare extremes heuristic
var1 = "relative_corr"
m = im1_p1900
var2 = "periods1900"
coef1 <- coef(m)[var1]
coef3 <- coef(m)[ paste0(var1,":", var2)]
med <- median(eval(parse(text = paste0("im1_p1900$model$", var2))))
hline <- coef1 + med*coef3

p1900_1 <- interplot(im1_p1900, var1 = "relative_corr", var2 = "periods1900", hist = T, point = T)+
  labs(y="Corruption", x="Year", title = 'Downes (2009) Win/Draw/Lose') +
  ggplot2::geom_hline(yintercept=hline, linetype="dashed") +
  ggplot2::geom_hline(yintercept=0) +
  theme_bw()


im2_p1900<- glm(wl ~ relative_corr*periods1900 +
                  initally + targally + v2x_libdem+ v2x_ex_military +concap + capasst+
                  qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                data = LEDs_data_f,  family = binomial(link = "probit"))
summary(im2_p1900)
im2sr_p1900<- coeftest(im2_p1900, 
                       vcov = sandwich(x = im2_p1900, adjust = TRUE),
                       cluster = war) #stata robust, slightly different because of of PC rather than HC

#grab relative_corr from interplot.medline for compare extremes heuristic
var1 = "relative_corr"
m = im1_p1900
var2 = "periods1900"
coef1 <- coef(m)[var1]
coef3 <- coef(m)[ paste0(var1,":", var2)]
med <- median(eval(parse(text = paste0("im1_p1900$model$", var2))))
hline <- coef1 + med*coef3

p1900_2 <- interplot(im1_p1900, var1 = "relative_corr", var2 = "periods1900", hist = T, point = T)+
  labs(y="Corruption", x="Year", title = 'Reiter and Stam (2002) Win/Lose') +
  ggplot2::geom_hline(yintercept=hline, linetype="dashed") +
  ggplot2::geom_hline(yintercept=0) +
  theme_bw()

im3_p1900<- glm(wldownes ~ relative_corr*periods1900 +
                  initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                  qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                data = LEDs_data_f,  family = binomial(link = "probit"))
summary(im3_p1900)

im3sr_p1900<- coeftest(im3_p1900, 
                       vcov = sandwich(x = im3_p1900, adjust = TRUE),
                       cluster = war) #stata robust, slightly different because of of PC rather than HC

#grab relative_corr from interplot.medline for compare extremes heuristic
var1 = "relative_corr"
m = im3_p1900
var2 = "periods1900"
coef1 <- coef(m)[var1]
coef3 <- coef(m)[ paste0(var1,":", var2)]
med <- median(eval(parse(text = paste0("im3_p1900$model$", var2))))
hline <- coef1 + med*coef3

p1900_3 <- interplot(im3_p1900, var1 = "relative_corr", var2 = "periods1900", hist = T, point = T)+
  labs(y="Corruption", x="Year", title = "Downes (2009) Win/Lose") +
  ggplot2::geom_hline(yintercept=hline, linetype="dashed") +
  ggplot2::geom_hline(yintercept=0) +
  theme_bw()


im4_p1900 <- lm(logLER ~ relative_corr*periods1900 +
                  + initally +targally + v2x_libdem+ v2x_ex_military +concap+ capasst+
                  qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
                data = LEDs_data_f) 
summary(im4_p1900)
im4sr_p1900<- coeftest(im4_p1900, 
                       vcov = vcovHC,
                       type = "HC1",
                       cluster = war) #stata robust 

#grab relative_corr from interplot.medline for compare extremes heuristic
var1 = "relative_corr"
m = im4_p1900
var2 = "periods1900"
coef1 <- coef(m)[var1]
coef3 <- coef(m)[ paste0(var1,":", var2)]
med <- median(eval(parse(text = paste0("im4_p1900$model$", var2))))
hline <- coef1 + med*coef3

p1900_4 <- interplot(im4_p1900, var1 = "relative_corr", var2 = "periods1900", hist = T, point = T)+
  labs(y="Corruption", x="Year", title = "Cochran and Long (2017) Battle Performance") +
  ggplot2::geom_hline(yintercept=hline, linetype="dashed") +
  ggplot2::geom_hline(yintercept=0) +
  theme_bw()


stargazer(im1_p1900, im2_p1900,im3_p1900, im4_p1900,
          se = list(im1sr_p1900[,2], im2sr_p1900[,2], im3sr_p1900[,2], im4sr_p1900[,2]),
          covariate.labels = c("Relative Corruption", "Post-1900", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4", "Corruption*Post-1900"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose', 
                            "Downes (2009) Win/Lose",
                            "Cochran and Long (2017) Battle Performance"
          ))

grid.arrange(p1900_1, p1900_2, p1900_3, p1900_4, nrow = 2, bottom = "Model")

#######Our own models 50 Year Periods Fixed Effects ########

im1_p50ff<- polr(as.factor(wdldownes) ~ relative_corr  + 
                   initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                   qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4 + as.factor(periods50), 
                 data = LEDs_data_f, method = c("probit"))
summary(im1_p50ff)
im1sr_p50ff<- coeftest(im1_p50ff, 
                       vcov = sandwich(x = im1_p50ff, adjust = TRUE),
                       cluster = war) #stata robust, slightly different because of of PC rather than HC

im2_p50ff<- glm(wl ~ relative_corr +
                  initally + targally + v2x_libdem+ v2x_ex_military +concap + capasst+
                  qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4 + as.factor(periods50), 
                data = LEDs_data_f,  family = binomial(link = "probit"))
summary(im2_p50ff)
im2sr_p50ff<- coeftest(im2_p50ff, 
                       vcov = sandwich(x = im2_p50ff, adjust = TRUE),
                       cluster = war) #stata robust, slightly different because of of PC rather than HC

im3_p50ff<- glm(wldownes ~ relative_corr +
                  initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                  qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4 + as.factor(periods50), 
                data = LEDs_data_f,  family = binomial(link = "probit"))
summary(im3_p50ff)

im3sr_p50ff<- coeftest(im3_p50ff, 
                       vcov = sandwich(x = im3_p50ff, adjust = TRUE),
                       cluster = war) #stata robust, slightly different because of of PC rather than HC


#with LEDs

im4_p50ff <- lm(logLER ~ relative_corr +
                  + initally +targally + v2x_libdem+ v2x_ex_military +concap+ capasst+
                  qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4 + as.factor(periods50), 
                data = LEDs_data_f) 
summary(im4_p50ff)
im4sr_p50ff<- coeftest(im4_p50ff, 
                       vcov = vcovHC,
                       type = "HC1",
                       cluster = war) #stata robust 


stargazer(im1_p50ff, im2_p50ff,im3_p50ff, im4_p50ff,
          se = list(im1sr_p50ff[,2], im2sr_p50ff[,2], im3sr_p50ff[,2], im4sr_p50ff[,2]),
          covariate.labels = c("Relative Corruption", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4", 
                               "1874-1923", "1924-1973","1974-2023"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose', 
                            "Downes (2009) Win/Lose",
                            "Cochran and Long (2017) Battle Performance"
          ))

#######Our own models Post-1900 Control ########

im1_p1900ff<- polr(as.factor(wdldownes) ~ relative_corr + 
                     initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
                     qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4 + as.factor(periods1900), 
                   data = LEDs_data_f, method = c("probit"))
summary(im1_p1900ff)
im1sr_p1900ff<- coeftest(im1_p1900ff, 
                         vcov = sandwich(x = im1_p1900ff, adjust = TRUE),
                         cluster = war) #stata robust, slightly different because of of PC rather than HC


im2_p1900ff<- glm(wl ~ relative_corr +
                    initally + targally + v2x_libdem+ v2x_ex_military +concap + capasst+
                    qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4 + as.factor(periods1900), 
                  data = LEDs_data_f,  family = binomial(link = "probit"))
summary(im2_p1900ff)
im2sr_p1900ff<- coeftest(im2_p1900ff, 
                         vcov = sandwich(x = im2_p1900ff, adjust = TRUE),
                         cluster = war) #stata robust, slightly different because of of PC rather than HC


im3_p1900ff<- glm(wldownes ~ relative_corr +
                    initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
                    qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4+ as.factor(periods1900), 
                  data = LEDs_data_f,  family = binomial(link = "probit"))
summary(im3_p1900ff)

im3sr_p1900ff<- coeftest(im3_p1900ff, 
                         vcov = sandwich(x = im3_p1900ff, adjust = TRUE),
                         cluster = war) #stata robust, slightly different because of of PC rather than HC


#with LEDs

im4_p1900ff <- lm(logLER ~ relative_corr +
                    + initally +targally + v2x_libdem+ v2x_ex_military +concap+ capasst+
                    qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4+ as.factor(periods1900), 
                  data = LEDs_data_f) 
summary(im4_p1900ff)
im4sr_p1900ff<- coeftest(im4_p1900ff, 
                         vcov = vcovHC,
                         type = "HC1",
                         cluster = war) #stata robust 



stargazer(im1_p1900ff, im2_p1900ff,im3_p1900ff, im4_p1900ff,
          se = list(im1sr_p1900ff[,2], im2sr_p1900ff[,2], im3sr_p1900ff[,2], im4sr_p1900ff[,2]),
          covariate.labels = c("Relative Corruption", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4", "Post-1900"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose', 
                            "Downes (2009) Win/Lose",
                            "Cochran and Long (2017) Battle Performance"
          ))


#######Model Fit Comparison#########
#Correlation Test between Vdems Democracy and Corruption Measure
cor.test(LEDs_data_f$v2x_libdem, LEDs_data_f$v2x_corr, method=c("pearson"))
# 
test <-LEDs_data_f %>% subset(select = c(wdldownes,wl,logLER,wldownes ))
test$wdldownes <- as.numeric(test$wdldownes)
cortest <- cor(test, use = "complete.obs")
stargazer(cortest)
#WL
LEDs_data_wl<- LEDs_data_f %>% drop_na(wl)

#Reiter and Stam (2002)  - no corruption 
iRSm1AIC <- AIC(iRSm1)
iRSm1BIC <- BIC(iRSm1)
pred_gam <-as.numeric(predict(iRSm1, type="response"))
p_gam <- ROCR::prediction(pred_gam, LEDs_data_wl$wl)
iRSm1pred<- performance( p_gam, "auc")@y.values[[1]]

#Reiter and Stam (2002)  - corruption 
iRSm2AIC <- AIC(iRSm2)
iRSm2BIC <- BIC(iRSm2)
pred_gam <-as.numeric(predict(iRSm2, type="response"))
p_gam <- ROCR::prediction(pred_gam, LEDs_data_wl$wl)
iRSm2pred<- performance( p_gam, "auc")@y.values[[1]]

#Own Model
im2AIC <- AIC(m2)
im2BIC <- BIC(m2)
pred_gam <-as.numeric(predict(m2, type="response"))
p_gam <- ROCR::prediction(pred_gam, LEDs_data_wl$wl)
im2pred <- performance( p_gam, "auc")@y.values[[1]]

#Reiter and Stam (2009)  - no corruption 
iRS2m1AIC <- AIC(iRS2m1)
iRS2m1BIC <- BIC(iRS2m1)
pred_gam <-as.numeric(predict(iRS2m1, type="response"))
p_gam <- ROCR::prediction(pred_gam, LEDs_data_wl$wl)
iRS2m1pred<- performance( p_gam, "auc")@y.values[[1]]

#Reiter and Stam (2009)  - corruption 
iRS2m2AIC <- AIC(iRS2m2)
iRS2m2BIC <- BIC(iRS2m2)
pred_gam <-as.numeric(predict(iRS2m2, type="response"))
p_gam <- ROCR::prediction(pred_gam, LEDs_data_wl$wl)
iRS2m2pred<- performance( p_gam, "auc")@y.values[[1]]

Model <- c("Proposed Model", "Reiter & Stam (2002)", "Reiter & Stam (2002)/corruption", "Reiter and Stam (2009)", "Reiter and Stam (2009)/corruption")
AIC <- c(im2AIC, iRSm1AIC, iRSm2AIC, iRS2m1AIC, iRS2m2AIC)
BIC <- c(im2BIC, iRSm1BIC, iRSm2BIC, iRS2m1BIC, iRS2m2BIC)
AUC <- c(im2pred, iRSm1pred, iRSm2pred, iRS2m1pred, iRS2m2pred)
wl_fit <- data.frame(Model, AIC, BIC, AUC)

#WDL
LEDs_data_wdl<- LEDs_data_f %>% drop_na(wdldownes)

#Downes2009 - no corruption 
iDm1AIC <- AIC(iDm1)
iDm1BIC <- BIC(iDm1)
pred_gam <-predict(iDm1, type="class")
pred_gam <- as.numeric(pred_gam)
Dfull <- preDATA(LEDs_data_wdl$wdldownes,pred_gam)
Dvec.full <- Dfull$Dvec
vus("full", T = pred_gam, Dvec = Dvec.full)
iDm1pred<- 0.3537
lipsitz.test(iDm1)
brant(iDm1)

#Downes2009 - corruption 
iDm2AIC <- AIC(iDm2)
iDm2BIC <- BIC(iDm2)
pred_gam <-predict(iDm2, type="class")
pred_gam <- as.numeric(pred_gam)
Dfull <- preDATA(LEDs_data_wdl$wdldownes,pred_gam)
Dvec.full <- Dfull$Dvec
vus("full", T = pred_gam, Dvec = Dvec.full)
iDm2pred <- 0.3582
lipsitz.test(iDm2)
brant(iDm2)


#Own Model
im1AIC <- AIC(m1)
im1BIC <- BIC(m1)
pred_gam <-predict(m1, type="class")
pred_gam <- as.numeric(pred_gam)
Dfull <- preDATA(as.factor(LEDs_data_wdl$wdldownes),pred_gam)
Dvec.full <- Dfull$Dvec
vus("full", T = pred_gam, Dvec = Dvec.full)
im1pred <-  0.3541
lipsitz.test(m1)
brant(m1)

Model <- c("Proposed Model", "Downes (2009)", "Downes (2009)/corruption")
AIC <- c(im1AIC, iDm1AIC, iDm2AIC)
BIC <- c(im1BIC, iDm1BIC, iDm2BIC)
VUS <- c(im1pred, iDm1pred, iDm2pred)
wdl_fit <- data.frame(Model, AIC, BIC, VUS)

#LED

#Own Model
im3AIC <- AIC(m3)
im3BIC <- BIC(m3)
im3AR2 <- summary(m3)$adj.r.squared

#Reiter and Stam (2002)  - no corruption 
iRSCLm1AIC <- AIC(iRSCLm1)
iRSCLm1BIC <- BIC(iRSCLm1)
iRSCLm1AR2 <- summary(iRSCLm1)$adj.r.squared

#Reiter and Stam (2002)  - corruption 
iRSCLm2AIC <- AIC(iRSCLm2)
iRSCLm2BIC <- BIC(iRSCLm2)
iRSCLm2AR2 <- summary(iRSCLm2)$adj.r.squared

#Downes2009 - no corruption 
iDCLm1AIC <- AIC(iDCLm1)
iDCLm1BIC <- BIC(iDCLm1)
iDCLm1AR2 <- summary(iDCLm1)$adj.r.squared

#Downes2009 - corruption 
iDCLm2AIC <- AIC(iDCLm2)
iDCLm2BIC <- BIC(iDCLm2)
iDCLm2AR2 <- summary(iDCLm2)$adj.r.squared

#Reiter and Stam (2009)  - no corruption 
iRSCL2m1AIC <- AIC(iRSCL2m1)
iRSCL2m1BIC <- BIC(iRSCL2m1)
iRSCL2m1AR2 <- summary(iRSCL2m1)$adj.r.squared

#Reiter and Stam (2009)  - corruption 
iRSCL2m2AIC <- AIC(iRSCL2m2)
iRSCL2m2BIC <- BIC(iRSCL2m2)
iRSCL2m2AR2 <- summary(iRSCL2m2)$adj.r.squared


Model <- c("Proposed Model", "Reiter & Stam (2002)", "Reiter & Stam (2002)/corruption",
           "Downes (2009)", "Downes (2009)/corruption", "Reiter & Stam (2009)", 
           "Reiter & Stam (2009)/corruption")
AIC <- c(im3AIC, iRSCLm1AIC, iRSCLm2AIC, iDCLm1AIC, iDCLm2AIC, iRSCL2m1AIC, iRSCL2m2AIC)
BIC <- c(im3BIC, iRSCLm1BIC, iRSCLm2BIC, iDCLm1BIC, iDCLm2BIC, iRSCL2m1BIC, iRSCL2m2BIC)
`Ajusted R^2`<- c(im3AR2, iRSCLm1AR2, iRSCLm2AR2, iDCLm1AR2, iDCLm2AR2, iRSCL2m1AR2, iRSCL2m2AR2)
led_fit <- data.frame(Model, AIC, BIC, `Ajusted R^2`)

#####Graphs

#wl_fit
wl_fit_AIC <- ggplot(data=wl_fit, aes(x= reorder(Model, -AIC) , y=AIC)) +
  geom_bar(stat = 'identity')+ 
  geom_text(aes(label=round(AIC, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal() + 
  labs(x="") 

wl_fit_BIC <-ggplot(data=wl_fit, aes(x= reorder(Model, -BIC) , y=BIC)) +
  geom_bar(stat = 'identity')+ 
  geom_text(aes(label=round(BIC, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+ 
  labs(x="") 

wl_fit_AUC <- ggplot(data=wl_fit, aes(x= reorder(Model, AUC) , y=AUC)) +
  geom_bar(stat = 'identity')+ 
  geom_text(aes(label=round(AUC, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+ 
  labs(x="") 


grid.arrange(wl_fit_AIC,wl_fit_BIC,wl_fit_AUC, nrow = 3, bottom = "Model")
#left = "Number of Groups", bottom = "Step")


#wdl_fit
wdl_fit_AIC <- ggplot(data=wdl_fit, aes(x= reorder(Model, -AIC) , y=AIC)) +
  geom_bar(stat = 'identity') + 
  geom_text(aes(label=round(AIC, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+ 
  labs(x="") 

wdl_fit_BIC <-ggplot(data=wdl_fit, aes(x= reorder(Model, -BIC) , y=BIC)) +
  geom_bar(stat = 'identity')+ 
  geom_text(aes(label=round(BIC, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+ 
  labs(x="") 

wdl_fit_VUS <-ggplot(data=wdl_fit, aes(x= reorder(Model, VUS) , y=VUS)) +
  geom_bar(stat = 'identity')+ 
  geom_text(aes(label=round(VUS, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+ 
  labs(x="") 

grid.arrange(wdl_fit_AIC,wdl_fit_BIC,wdl_fit_VUS, nrow = 3, bottom = "Model")
#left = "Number of Groups", bottom = "Step")

#led_fit
led_fit_AIC <- ggplot(data=led_fit, aes(x= reorder(Model, -AIC) , y=AIC)) +
  geom_bar(stat = 'identity')+ 
  geom_text(aes(label=round(AIC, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+ 
  labs(x="") 

led_fit_BIC <- ggplot(data=led_fit, aes(x= reorder(Model, -BIC) , y=BIC)) +
  geom_bar(stat = 'identity')+ 
  geom_text(aes(label=round(BIC, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+ 
  labs(x="") 

led_fit_AR2 <-ggplot(data=led_fit, aes(x= reorder(Model, Ajusted.R.2) , y=Ajusted.R.2)) +
  geom_bar(stat = 'identity')+ 
  geom_text(aes(label=round(Ajusted.R.2, digits = 3)), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+ 
  labs(x="") 

grid.arrange(led_fit_AIC,led_fit_BIC,led_fit_AR2, nrow = 3, bottom = "Model")
#left = "Number of Groups", bottom = "Step")

